%% dXt=int SIGMAt W sigmat=W or OU
%the revisible part are K,N,RNK, R(s), sigma(s), y0, tspan
clear all
casenum=5;
% index J for X
alphaK=20; alphas=[zeros(1,alphaK);eye(alphaK);2*eye(alphaK);index2(alphaK)];
fname='alphas'; save(fname,'alphas');
% index J for sigma
betas = [zeros(1,alphaK);eye(alphaK)];%...
fname='betas';  save(fname,'betas');
% time T
T = 1/2;
% construct orthogoal basis in L2[0,t]
orthbasis(T);
% initial condition for sigma
sig0=zeros(size(betas,1),1); sig0(1)=1;
% time mesh 
num=50;
tspan = linspace(0,T,num+1);
% ms' value
msvalue(num);
load('msval');
% solver of sigma
options = odeset('AbsTol',1e-12,'RelTol',1e-6);
% % % %parameters of ou
a=1;m=1;b=1; para=[a,m,b,num];
disp('solving sigma')
tic
%[time,sig] = ode45(@randomOU,tspan,sig0,options,para);
[time,sig] = ode45(@randomOU,tspan,sig0,[],para);
toc
% sig=zeros(num+1,betaK+1);
% sig(:,1)=1;
fname='sigma';
save(fname,'sig');
row=size(alphas,1);
X0=zeros(row,1);
% dt=0.5/num;
% solver of X
disp('solving X')
tic
%[time,X] = ode45(@randomtest,tspan,X0,options,[num,alphaK],casenum);
[time,X] = ode45(@randomtest,tspan,X0,[],[num,alphaK],casenum);
toc
% wsq{1}=@(s) s.^3/3/t;
% for k=2:alphaK
%     wsq{ii}=@(s) 2*t/(k-1)^2/pi^2*sin((k-1)*pi/t*s);
% end
% for tt=2:num+1
%     tmp=sum(sig(:,2:alphaK+1).*wsval(:,1:alphaK),2);
%     tmp=@(s) wsq{1}(s)
%     for jj=1:alphaK
%         tmp=@(s)
%     X(tt,1)=midptrule(tmp(1:tt),dt);
% end
% for ii=2:alphaK+1
%     for tt=2:num+1
%     tmp=sig(:,1).*wsval(:,ii-1);
%     X(tt,ii)=quad()%midptrule(tmp(1:tt),dt);
%     end
% end
% begins=alphaK+2; ends=2*alphaK+1;
% for ii=begins:ends
%     tmp=sqrt(2)*sig(:,ii-alphaK).*wsval(:,ii-alphaK-1);
%     X(tt,ii)=midptrule(tmp(1:tt),dt);
% end
% %X(:,alphaK+2:2*alphaK+1)=sqrt(2)*sig(:,2:alphaK+1).*wsval(:,1:alphaK);
% begins=2*alphaK+2; ends=begins+alphaK*(alphaK-1)/2-1;
% for ii=begins:ends
%     myk=find(alphas(ii,:)==1);
%     tmp=wsval(:,myk(1)).*sig(:,myk(2)+1)+wsval(:,myk(2)).*sig(:,myk(1)+1);
%     X(tt,ii)=midptrule(tmp(1:tt),dt);
% end

varitest(time',X,a,m,b,sig0(1),X.*X,casenum)



%% dXt=-a*m* int W dt
%the revisible part are K,N,RNK, R(s), sigma(s), y0, tspan
clear all
casenum=4;
% index J for X
alphaK=50; alphas=[zeros(1,alphaK);eye(alphaK)];
fname='alphas'; save(fname,'alphas');
% time T
T = 1/2;
% construct orthogoal basis in L2[0,t]
orthbasis(T);
% time mesh 
num=50;
% ms' value
msvalue(num);
load('msval');
%
a=1;m=1;b=1;
row=size(alphas,1);
X=zeros(num+1,row);
disp('solving X')
tic
for ii=2:alphaK+1
    X(:,ii)=-a*m*nsval(:,ii-1);
end
toc
time=linspace(0,0.5,num+1);
varitest(time',X,a,m,b,1,X.*X,casenum)


%% dXt=SIGMAt W sigma=OU or W
%the revisible part are K,N,RNK, R(s), sigma(s), y0, tspan
clear all
casenum=3;
% index J for X
alphaK=20; alphas=[zeros(1,alphaK);eye(alphaK);2*eye(alphaK);index2(alphaK)];
fname='alphas'; save(fname,'alphas');
% index J for sigma
betas = [zeros(1,alphaK);eye(alphaK)];%...
fname='betas';  save(fname,'betas');
% time T
T = 1/2;
% construct orthogoal basis in L2[0,t]
orthbasis(T);
% initial condition for sigma
sig0=zeros(size(betas,1),1); sig0(1)=0;
% time mesh 
num=50;
tspan = linspace(0,T,num+1);
% ms' value
msvalue(num);
load('msval');
% solver of sigma
options = odeset('AbsTol',1e-12,'RelTol',1e-6);
% % % %parameters of ou
a=0;m=1;b=1; para=[a,m,b,num];
disp('solving sigma')
tic
%[time,sig] = ode45(@randomOU,tspan,sig0,options,para);
[time,sig] = ode45(@randomOU,tspan,sig0,[],para);
toc
% sig=zeros(num+1,betaK+1);
% sig(:,1)=1;
fname='sigma';
save(fname,'sig');
row=size(alphas,1);
X=zeros(num+1,row);
disp('solving X')
tic
X(:,1)=sum(sig(:,2:alphaK+1).*wsval(:,1:alphaK),2);
for ii=2:alphaK+1
    X(:,ii)=sig(:,1).*wsval(:,ii-1);
end
X(:,alphaK+2:2*alphaK+1)=sqrt(2)*sig(:,2:alphaK+1).*wsval(:,1:alphaK);
begins=2*alphaK+2; ends=begins+alphaK*(alphaK-1)/2-1;
for ii=begins:ends
    myk=find(alphas(ii,:)==1);
    X(:,ii)=wsval(:,myk(1)).*sig(:,myk(2)+1)+wsval(:,myk(2)).*sig(:,myk(1)+1);
end
toc
varitest(time',X,a,m,b,sig0(1),X.*X,casenum)


%% dXt=SIGMAtdW2
%the revisible part are K,N,RNK, R(s), sigma(s), y0, tspan
clear all
casenum=2;
% index J for X
alphaK=50; alphas=[zeros(1,alphaK);eye(alphaK);2*eye(alphaK);index2(alphaK)];
fname='alphas'; save(fname,'alphas');
% index J for sigma
betaK=50;  betas = [zeros(1,betaK);eye(betaK)];%...
fname='betas';  save(fname,'betas');
% time T
T = 1/2;
% construct orthogoal basis in L2[0,t]
orthbasis(T);
% initial condition for sigma
sig0=zeros(size(betas,1),1); sig0(1)=0;
% time mesh 
num=50;
tspan = linspace(0,T,num+1);
% ms' value
msvalue(num)
% solver of sigma
options = odeset('AbsTol',1e-12,'RelTol',1e-6);
% % % %parameters of ou
a=0;m=0;b=1; para=[a,m,b,num];
disp('solving sigma')
tic
[time,sig] = ode45(@randomOU,tspan,sig0,[],para);
toc
% sig=zeros(num+1,betaK+1);
% sig(:,1)=1;
fname='sigma';
save(fname,'sig');
%calculate value table
load('msval')
dt=1/100;
intwi=zeros(alphaK,num+1); intwisig0=zeros(alphaK,num+1);
intwisigj=zeros(num+1,alphaK,betaK);
disp('numerical integration')
tic
for ii=1:alphaK
    for tt=1:num+1
        intwi(ii,tt)=midptrule(wsval(1:tt,ii),dt);
        intwisig0(ii,tt)=midptrule(wsval(1:tt,ii).*sig(1:tt,1),dt);
    end
end
for ii=1:alphaK
    for jj=1:betaK
        for tt=1:num+1
            intwisigj(tt,ii,jj)=midptrule(wsval(1:tt,ii).*sig(1:tt,1+jj),dt);
        end
    end
end
toc
% solver of X
betarow=size(betas,1);
row=betarow*(alphaK+1)+alphaK+alphaK*(alphaK-1)/2;
X=zeros(num+1,row);
X(:,1)=-1/2*b*time';
for ii=2:alphaK+1
    begins=(ii-1)*(betaK+1)+1;
    X(:,begins)=-a*m*intwi(ii-1,:)'+a*intwisig0(ii-1,:)';
    for jj=1:betaK
        X(:,begins+jj)=sig(:,jj+1).*wsval(:,ii-1)+a*intwisigj(:,ii-1,jj);
    end
end
begins=(alphaK+1)*(betaK+1)+1; ends=(alphaK+1)*(betaK+2);
% X(:,:)=0;
% X(:,1)=time;
for ii=begins:ends
    X(:,ii)=-1/2*b*sqrt(2)*wsval(:,ii-begins+1).^2;
end
begins=ends+1; ends=ends+(alphaK-1)*alphaK/2; shift=(alphaK+1)*betaK+1;
for ii=begins:ends
    myk=find(alphas(ii-shift,:)==1);
    X(:,ii)=-b*wsval(:,myk(1)).*wsval(:,myk(2));
end


disp('sloving X')

% error
%meanerr(time',y,R,sig,y0,Y/numpath)
varitest(time',X,a,m,b,sig0(1),X.*X,casenum)
%varXOU(time',X,a,m,b,X0(1),X.*X)
%% dXt=WdW Yt=2(Xt+t/2)
casenum=1;
% index J for Yt
betaK=50;  betas = [zeros(1,betaK);eye(betaK);2*eye(betaK);index2(betaK)];
% initial condition for Yt
betarow=size(betas,1);
row=betarow;
X=zeros(num+1,row);
% time mesh 
num=50;
tspan = linspace(0,T,num+1)';
X(:,1)=tspan;
T = 1/2;
% construct orthogoal basis in L2[0,t]
orthbasis(T);
% calculate value table
msvalue(num);
% loading
load('msval');
disp('solving X')
tic
begins=betaK+2; ends=2*betaK+1;
for ii=begins:ends
    X(:,ii)=sqrt(2)*wsval(:,ii-betaK-1).^2;
end
begins=2*betaK+2; ends=2*betaK+1+betaK*(betaK-1)/2;
for ii=begins:ends
    myk=find(betas(ii,:)==1);
    X(:,ii)=2*wsval(:,myk(1)).*wsval(:,myk(2));
end
toc
varitest(time',X,1,1,1,1,X.*X,casenum)